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Generally "exact" Quantum Monte Carlo computations for the ground state of many Bosons make 
use of importance sampling. The importance sampling is based, either on a guiding function or on 
an initial variational wave function. Here we investigate the need of importance sampling in the case 
of Path Integral Ground State (PIGS) Monte Carlo. PIGS is based on a discrete imaginary time 
evolution of an initial wave function with a non zero overlap with the ground state, that gives rise 
to a discrete path which is sampled via a Metropolis like algorithm. In principle the exact ground 
state is reached in the limit of an infinite imaginary time evolution, but actual computations are 
based on finite time evolutions and the question is whether such computations give unbiased exact 
results. We have studied bulk liquid and solid 4 He with PIGS by considering as initial wave function 
a constant, i.e. the ground state of an ideal Bose gas. This implies that the evolution toward the 
ground state is driven only by the imaginary time propagator, i.e. there is no importance sampling. 
For both the phases we obtain results converging to those obtained by considering the best available 
variational wave function (the Shadow wave function) as initial wave function. Moreover we obtain 
the same results even by considering wave functions with the wrong correlations, for instance a 
wave function of a strongly localized Einstein crystal for the liquid phase. This convergence is 
true not only for diagonal properties such as the energy, the radial distribution function and the 
static structure factor, but also for off-diagonal ones, such as the one-body density matrix. This 
robustness of PIGS can be traced back to the fact that the chosen initial wave function acts only at 
the beginning of the path without affecting the imaginary time propagator. From this analysis we 
conclude that zero temperature PIGS calculations can be as unbiased as those of finite temperature 
Path Integral Monte Carlo. On the other hand, a judicious choice of the initial wave function greatly 
improves the rate of convergence to the exact results. 



I. INTRODUCTION 

Among the available methods to investigate the prop- 
erties of strongly interacting many-body quantum sys- 
tems, Quantum Monte Carlo (QMC) ones hold a rele- 
vant position since some of the QMC methods can pro- 
vide "exact" expectation values. As for all the Monte 
Carlo methods, QMC results are affected by statistical 
uncertainties, but the property that makes some of them 
"exact" is the possibility of reducing within this unavoid- 
able errors all the systematic errors introduced by the in- 
volved approximations. This is true for Boson system at 
zero and at finite temperature, but the studies of Fcrmion 
systems, excited states or real time dynamics all suffer 
from sign problems which still have precluded the devel- 
opment of such kind of "exact" methods. A first great 
subdivision of "exact" QMC methods is between finite 
and zero temperature methods. Among T — K meth- 
ods we find Green Function Monte Carlo (GFMC), 1 ^ 
Diffusion Monte Carlo (DMC), 4 Path Integral Ground 
State Monte Carlo (PIGS)£ and Reptation Monte Carlo 
(RMC)j2 At finite temperature, the major role is played 
by Path Integral Monte Carlo (PIMC)^ 

The main difference between finite and zero tempera- 
ture methods is that ground state methods rely, to dif- 
ferent extents, on a trial wave function for the impor- 
tance sampling, while in PIMC no such wave function 
is involved at all. PIMC needs only information on the 
interaction among particles as input, and the main dif- 
ficulty to overcome, beyond the propagator accuracy, is 
the sampling ergodicity. Recently a great step forward 



in this direction has been realized with the advent of the 
worm algorithm^ The first GFMC computation did not 
use importance sampling^ but this computation was for 
a small system of 32 particles. All the other computa- 
tion at zero temperature have used a model wave function 
that was given as input. With respect to PIGS method, if 
, in principle, convergence can be achieved for any initial 
wave function that has a finite overlap with the ground 
state, the question is if convergence can be achieved in 
a real computation. Some progress have been achieved 
recently both for a finite 9 and for a bulk Boson system.— 
For example, a Jastrow wave function, that describes the 
bulk phase of a quantum liquid, was employed as ini- 
tial wave function in the study of small parahydrogen 
clusters^ In the case of a two dimensional crystal of 4 Hc, 
convergence to the same result was found in the com- 
putation of the one-body density matrix starting from 
two initial wave functions with "opposite" properties^ 
one with Bose-Einstein condensate and the other with- 
out. Stimulated by these results, we have undergone a 
systematic check for a realistic Hamiltonian for 4 He and 
we have found that PIGS methods converge to the ex- 
act result regardless of the chosen initial wave function. 
Even better, PIGS converges even if the wave function 
contains wrong correlations or no correlations at all. 

As a test-bed we have considered the bulk liquid and 
solid phases of a strongly interacting boson system such 
as 4 He. By projecting very different wave functions, we 
find converging results for diagonal properties like the en- 
ergy, the radial distribution function, the static structure 
factor and also for off-diagonal properties like the one- 
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body density matrix. For both the liquid and the solid 
phases we consider a shadow wave function (SWF)jii 
which is the best available variational wave function 12 
for 4 He systems, and the constant wave function (CWF), 
which is the ground state wave function of the ideal Bosc 
gas and does not contain any correlation. The latter cor- 
responds to have no importance sampling since the initial 
wave function gives the same weight to any configuration 
of the particles. Furthermore, for the liquid phase, we 
have considered also a localized wave function which de- 
scribes an Einstein solid. Our question is whether PIGS 
is able to recover the exact ground state starting from 
really different systems, i.e. an ideal gas or an Einstein 
crystal for the quantum liquid and an ideal gas for the 
quantum solid. We find that for all the considered ini- 
tial wave functions PIGS converges to the exact ground 
state. Only the convergence rate depends on the initial 
wave function, and turns out to be slower for the worse 
ones. 

The paper is organized as follows: Scc.|ll]dcals with the 
description of the PIGS method, of the test bed systems 
and of the initial wave functions. Results are presented 
and discussed in Sec. IIIIl Sec. IIVI contains our conclu- 
sions. 



II. METHOD 

A. The PIGS method 

It is well known that, given a Hamiltonian H for a 
quantum system of N particles, the ground state wave 
function if}$ in the position representation can be ob- 
tained as the t — > oo limit of an imaginary time (t = 
it/h) evolution of an initial (trial) wave function ipT, pro- 
vided that (i/jq\iPt) 7^ 0: 



ipa{R) — lim 



e - T (ff-E )^ T 



(1) 



Regardless of a normalization constant, which is not in- 
volved in the Monte Carlo sampling and thus will be 
dropped in the following, an accurate approximation of 
the ground state wave function is given by 



if) T ( R ) = J dR ' G(R,R';T)ip T (R). 



(2) 



This originates from the action of the imaginary time 
projector G = e~ rH , which exponentially removes from 
ipT any overlap with the excited states during the imag- 
inary time evolution. The evolved wave function ip T and 
Eq. (fT]) provide the basis for all the zero temperature 
QMC methods. 

A first problem rises: the imaginary time propagator 
G can be accurately written only for small values of t, 
but a large r limit is necessary to ensure the convergence 
to ipQ. One possible strategy to overcome this problem is 



to reach the large r limit by means of a recursive proce- 
dure, as for example in GFMCi and DMCi^ Both these 
methods reach an extremely accurate approximation of 
the ground state wave function multiplied by -0T by it- 
erating equation (fT]) by means of random walks. With 
these methods the trial wave function has a strategical 
role since it is involved at each iteration step, where it is 
used as importance sampling function. 

On the contrary, the PIGS method^ interprets the 
imaginary time propagator as a density matrix opera- 
tor corresponding to an inverse temperature /3 = r, and 
then, by exploiting the factorization property 

G(R,R';t = t 1 +t 2 ) = J dR" G(R,R";t 1 )G{R",R';t 2 ) 

(3) 

the large r propagator G is written as a convolution of 
small imaginary time propagators G(R, R'; St) as in stan- 
dard path integral formalism.— Thus the density matrix 
operator is broken up into M small pieces with a time 
step St = t/M and the approximated ground state wave 
function reads 

4> T {R)= J dR 1 ---dR M G{R,R 1 ;ST)G{R 1 ,R 2 ;ST)--- 



x G(R M -i,Rm;St)iPt(Rm)- 



(4) 



An appealing feature, peculiar of this method is that, in 
ij} T , the ansatz on ipT acts only at the starting point be- 
ing the full imaginary time path governed by G, which 
depends only on H. Once fixed St, the elementary evo- 
lution step in imaginary time is obtained acting with 
G(R,R';5t) on a quantum state; this action is usually 
called projection step due to the resulting increased over- 
lap with the ground state. The PIGS method reaches 
convergence to the ground state when adding further 
imaginary time projections to ip r does not provide any 
appreciable change of the results, i.e. the expectation 
values computed with ip T+n g T is compatible with those 
computed with ip T within the statistical errors for any 
integer n > 0. Since the convergence on the ground state 
for a finite system is exponentially fast^ the number M 
of required imaginary time projections is usually limited. 
Nevertheless M depends on the specific expectation value 
one is computing, thus a separate analysis of the conver- 
gence should be carried on for every computed quantity. 
The convergence rate is determined essentially by the 
quality of the initial wave function. The total imaginary 
time r required to clean the initial wave function up from 
excited state contributions is smaller for good (i.e. with 
large overlap with the ground state) initial wave func- 
tions. Also the accuracy of the approximation used for 
the small time G has a role in determining the computa- 
tional cost of the projection procedure: in fact, with bet- 
ter approximations of the small imaginary time G, one 
needs a smaller number of projection steps to reach con- 
vergence since larger time step St are allowed. Moreover 
the imaginary time projection procedure should be able, 
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in principle, to incorporate in the final wave function all 
the correct correlations, even if these were absent in ipx- 
Thus it should be possible, in principle, to obtain the 
correct ground state wave function, if r is large enough, 
even if ipT is extremely poor, and the results do not suf- 
fer from any variational bias. The intent of the present 
study is to show that this is indeed the case with values 
of t manageable with present computational resources. 

Like for finite temperature PIMC^ also within PIGS 
the quantum system is mapped into a system of classical 
polymers, but in this case the polymers start and end 
on the initial wave function ipx ^ ' 14 i 15 resulting in open 
linear polymers instead of the close ring ones of PIMC. 7 
Each open polymer represents the full imaginary time 
path of a quantum particle that is sampled by means 
of the Metropolis algorithm. Thus, the entire imaginary 
time evolution of the system is sampled at each Monte 
Carlo step, contrary to what happens for DMC, for ex- 
ample. Expectation values of an operator O reads 

6 _ ^t\G(t)6G{t)\^ t ) 
h ^ t \G(t)G(t)\^ t ) 

and, if r is large enough to allow convergence to the 
ground state, (0) T is the ground state expectation value, 
without needing any extrapolation. However an analy- 
sis of the convergence as a function of r is needed, as 
explained before. 

The ground state expectation value of the energy can 
be obtained in several ways: the most largely employed 
one is the mixed estimate 

^ = (^t\G(2t)H\^ t ) + (^t\HG{2t)\^ t ) 
2(^ t |G(2t)|^t) 

Notice that in PIGS this mixed estimate is exact. In 
fact, since the Hamiltonian operator H commutes with 
the imaginary time evolution operator G, it is possible 
to obtain unbiased expectation value of the Hamiltonian 
inserting H at one of the ends of the path. If ipx is an 
accurate wave function, this estimator is preferable with 
respect to the direct one given by Eq. ([5|) with O = H, 
because it has typically a lower variance. This is no more 
true, however, when the initial wave function is particu- 
larly poor. In this case the fluctuations in the expecta- 
tion value of the Hamiltonian become so sizable that © 
is unusable. The direct estimator is seldom used because 
one has to compute the first and the second derivative 
of G(R,R';St); but this is not a serious problem if the 
imaginary time propagator is available in an analytical 
formulation, like the one we employed here. When ipT is 
a good wave function, as for the SWF case, we have veri- 
fied that mixed and direct estimators provide compatible 
results, and we report here the mixed one for the above 
mentioned reason. On the contrary, for poor wave func- 
tions, like CWF, we had to consider the direct estimator. 

As far as 4 He systems are concerned, because of the 
Bose statistic obeyed by the atoms, one has, in princi- 
ple, to account for permutations in the propagator G^^ 



Permutation moves are not strictly requested whenever 
the initial wave function has the correct Bose symme- 
try. In fact, the polymer configuration resulting after a 
permutation can be in principle reached with a combina- 
tion of standard sampling moves, since the polymers are 
openi 14 ! 15 This is not the case if i\>t is not Bose symmet- 
ric like the Gaussian wave function: permutation cycles 
among particles must be introduced in the sampling in 
order to get the exact ground state. For all xjjx we have 
implemented permutation sampling following Ref. [lrl In 
fact, implementing the sampling of permutations even 
for a Bose symmetric ipr turns out to greatly improve 
the ergodicity of the sampling. In addition, we have 
used swap moves^ because they increase the sampling 
efficiency when computing off-diagonal properties^ 



B. Test systems 

In order to test the PIGS method convergence proper- 
ties we have considered two bulk phases of a many-body 
strongly interacting Boson system: liquid and solid 4 He. 
Dealing with low temperature properties, 4 He atoms are 
described as structureless zero-spin bosons, interacting 
through a realistic two-body potential, that we assume 
to be the HFDHE2 Aziz potential^ 7 -. For the liquid phase, 
we have considered a cubic box with periodic boundary 
conditions, containing N = 64 atoms at the equilibrium 
density pi = 0.0218A -3 . For the solid phase we have 
considered a cubic box with periodic boundary condi- 
tions designed to house a fee crystal of N = 32 atoms 
at the density p s = 0.0313A" 3 . In both cases we add 
standard tail corrections to the potential energy to ac- 
count for the finite size of the system by assuming the 
medium homogeneous (i.e. g(r) — 1) beyond L/2, where 
L is the size of the box. Obviously, this is not an ac- 
curate assumption specially for the solid phase in such 
a small box, but our main purpose here is to show that 
PIGS method is able to reach the same results indepen- 
dently on the considered initial wave function. Moreover 
we have studied the fee lattice, which is stabilized by the 
cell geometry and by the periodic boundary conditions, 
whereas at zero temperature, the hep lattice is the stable 
lattice. Computations of ground state properties of bulk 
4 He with accurate tail corrections can be found in the 
current literature i 12 i 18 



C. Initial wave functions 

The standard initial wave functions commonly used 
within PIGS method 5 - are the variational Jastrow wave 
function (JWF) for the liquid and the Jastrow-Nosanow 
(J-NWF) for the solid. A JWF represents the simplest 
possible choice of wave function for strongly interact- 
ing Bosons 1 ^ and it contains only two-body correlations. 
Using a McMillan pseudopotential, 20 the unnormalized 
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JWF reads as 

^j WF (i?)= ]J e~Hf«; . (7) 

Kj = l 

The physical meaning of this JWF is that, due to the 
sharp repulsive part of the interaction potential V in the 
Hamiltonian H , He atoms prefer to avoid each other. In 
the J-NWF the JWF is multiplied by a term like the one 
in Eq. (fT0|) below, that localizes the particles in a crys- 
talline order. In this work, however, in order to explore 
the convergence properties of the PIGS method, we have 
considered two wave functions of "opposite" quality: the 
best available one, that is the shadow wave function, and 
the poorest imaginable one, i.e. the constant wave func- 
tion. JWF will be considered only when computing the 
one-body density matrix in the liquid phase. 

The constant wave function is the ground state wave 
function of the ideal Bose gas, 

V-cwfOR) = 1- (8) 

It carries no correlation at all. We choose this wave func- 
tion because, allowing an unrestricted sampling of the full 
configurational space, it results in no importance sam- 
pling. Then the whole imaginary time projection proce- 
dure is driven only by the short time evolution operator 
G, without any input, and then any bias, from the initial 
state. Thus at the starting point the system is made up 
by free particles; if after a long enough imaginary time 
projection, PIGS turns out to be able to reach a strong 
correlated quantum liquid and quantum crystal by itself 
we can safely believe that no variational bias affects PIGS 
results. 

On the other hand, we choose as tpx a SWF given 
by variational computation in order to have as refer- 
ence results the ones coming from the projection of an 
initial wave function that is more accurate as possible, 
i.e. from a wave function whose overlap with the ex- 
act ground state is known to be large. In the SWF, 
additional correlations besides the standard two body 
terms are introduced via auxiliary variables which are 
integrated outM- This is done so efficiently that the crys- 
talline phase emerges as a spontaneously broken symme- 
try process, induced by the inter-particles correlations 
as the density is increased, without the need of any a 
priori knowledge of the equilibrium positions and with- 
out losing the translationally invariant form of the wave 
function. Thus SWF is able to describe both the liq- 
uid and the solid phase with the same functional form 
and it is explicitly Bose symmetric. The standard SWF 
functional form reads 

^swf(R) = MR) [ dSK(R,S)<t> s (S) (9) 

where S = (s\, s 2 , . . . , spj) is the set of auxiliary shadow 
variables, 4> r (R) is the standard Jastrow two body cor- 
relation term Q, K(R,S) is a kernel coupling each 



shadow to the corresponding real variable, and i/) s (S) is 
another Jastrow term describing the inter-shadow cor- 
relations. As usual,— we take K(R, S) Gaussian and 
in <fis(S) we use the rescaled and dilated He-He poten- 
tial V as pseudopotential. The variational parameters 
we use were chosen in order to minimize the expecta- 
tion value of the Hamiltonian H and are reported in 
Ref. [U Nowadays the SWF represents the best avail- 
able variational wave function for 4 He systems^ Re- 
cently, we have estimated-^ that, when describing a two 
dimensional solid, SWF overlap per particle with the 
true ground state is of about 99.8%, which ensures a 
fast convergence rate when projected within the PIGS 
method. The properties of the SWF are so peculiar that 
the PIGS method that has a SWF as ipT deserves an 
its own name and is dubbed SPIGS: Shadow Path Inte- 
gral Ground State method . 14 ' 15 In the picture of linear 
polymers, the presence of the shadow variables adds two 
extra variational links, one at each end of the polymer. 

In order to test how robust PIGS is, we consider also 
a "wrong" wave function: for the liquid phase we con- 
sider a Gaussian wave function, where each particle is 
harmonically localized around fixed positions {r*oi} 

N 

^cwF^n 6 ^ 1 ™' 2 ' ( io ) 

1=1 

i.e. ipT it the wave function of an Einstein harmonic solid. 
The parameter C = 8 A -2 is arbitrary and it is was cho- 
sen to ensure a strong localization of the particles around 
the positions {foi} that were taken over a regular cubic 
lattice within the simulation box. This wave function 
is evidently not translationally invariant and not Bose 
symmetric. Furthermore it does not contain any correla- 
tion between the particles, and all the information that 
it carries is that of a crystalline system, i.e. GWF is an 
extremely poor wave function for the liquid phase. This 
wrong initial wave function will provide a stringent test 
on the convergence properties of the PIGS methods. 

As far as the one-body density matrix computation in 
the liquid phase is concerned, the values of the parame- 
ters b and m in the JWF have been chosen equal to the 
ones of the corresponding Jastrow term in the SWF. 

D. Small time propagator 

One of the fundamental elements of path integral pro- 
jection Monte Carlo methods is the imaginary time prop- 
agator G, whose accuracy turns out to be crucial to the 
convergence to the exact results. The functional form 
of G for a generic r is unfortunately not known with 
exception of few particular cases, such as, for example, 
the free particle and the harmonic oscillator, but accu- 
rate approximations of G are obtainable in the small 
t regime /' 13 ' 22 In this work, we have chosen the Pair- 
Suzuki approximation 23 for the imaginary time propaga- 
tor, which is a pair-approximation of the fourth-order 
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Suzuki-Chin density matrix^ 

The Suzuki-Chin approximation is based on the fol- 
lowing factorization of the density matrix: 



-2StH 



e -^V ee -5rf e -^V Ce -Srf e -^V e 



(11) 



where T is the kinetic operator and V e and V c are given 
by 



V e = V + 



X 2 \ N 

6 *=i 



and 



V C = V + 



(1 - a)Sr 2 X 
6 



N 



EC*) 



(12) 



(13) 
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respectively, with V the potential operator, a an ar- 
bitrary constant in the range [0,1], A = h 2 /2m and 
Pi = V \V . The resulting imaginary time propagator 
is accurate to order (5r 4 , and has been successfully ap- 
plied to liquid He in two and three dimensions^ This 
approximation offers also the advantage that adjusting 
the parameter a it is possible to optimize the conver- 
gence, and a standard choice for a quantum system is 
a = 0.— A strategy to obtain a simpler, but equally accu- 
rate, approximation consists in applying a pair product 
assumption. 2 - For sufficiently short time steps, in fact, 
the many-body propagator (in imaginary time) is well 
approximated by the product of two-body propagators. 7 
In this approximation, the small time propagator reads 



G(R m: Rm +1 ;5T) = (4tt\6t 

N 

ri exp 



-3JV/2 



i,m Ti,m-\-l j 

AXSt 



exp (—u(nj jm , rjj.m+i)) 



(14) 



where u is given as 



u{r m ,r m+ i) 



4r [v e {r m ) + 2u c (r m+ i)] m odd 



St 



3 [2v c (r m ) + v e (r m+1 )] m even. 



The potentials v e (r) and v c (r) arc defined as 

2 



v e (r) = V(r) + cx-St 2 X 



dV 
Or 



2 {r) = V{r) + {l-a)-8T 2 \ 



&V 
Or 



(15) 



(16) 



where V(r) is the potential experienced by two 4 He 
atoms at a distance r. The advantage is that there is 
no need to calculate Fi. As for the full Suzuki-Chin 
approximation,— also for the Pair-Suzuki the operators 
corresponding to physical observables must be inserted 
only on odd time slices in the imaginary time path. 



Figure 1: Energy per 4 He atom E(t) vs. imaginary time 
step St. The total projection time is r = 0.1 K _1 . The 
calculations were carried out by projecting a SWF and a CWF 
for a system of 64 particles at the equilibrium density p = 
0.0218 A' 3 . Dashed lines are quartic fits to the data. Error 
bars, when not shown, are smaller than the used symbols. 



In order to fix the optimal small imaginary time step 
value, we have performed PIGS simulations with differ- 
ent initial wave functions. By considering decreasing St 
values with a fixed total projection time, r, we have taken 
the energy per particle E(r) as observable of reference. 
As an example, our results for SWF and CWF in the 
liquid phase are plotted in Fig. [T] We choose as optimal 
value St — 1/640 K _1 ; in fact, further reductions do not 
change the energy in a detectable way, i.e. within the 
statistical uncertainty. In Fig. Q] SWF and CWF do not 
converge to the same value simply because the consid- 
ered total projection time r in this test is not enough to 
ensure convergence of E(t) to the ground state energy 
for CWF (see Fig. Similarly, in the solid phase we 
take St — 1/960 K _1 . 



III. RESULTS 

Once set the optimal St value, we have computed the 
diagonal properties of the system for increasing total pro- 
jection time r until we reached convergence to a value 
that corresponds to the exact ground state result both 
for the liquid and for the solid phase. In the liquid phase 
we have computed also the one-body density matrix. 



A. Liquid 

1. PIGS results without importance sampling 

For the liquid phase we have projected a SWF and a 
CWF. The energy per particle as a function of the total 
projection time r for both the wave functions is plotted 
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Figure 2: Energy per particle £ as a function of the total 
projection time r obtained from PIGS simulations for liquid 
4 He at the equilibrium density p = 0.0218 A -3 by projecting 
a SWF (filled circles) and a CWF (open circles) and a GWF 
(open diamonds), r = result (filled circle) corresponds to 
the SWF variational estimate of E, the r = for the GWF 
is E = 122.08 ± 0.06 K and for CWF E is essentially infinite. 
Error bars are smaller than the used symbols. Dotted line 
indicates the convergence value E = —7.17 ± 0.02 K. 



in Fig. O We find that the energy converges, indepen- 
dently from the considered initial wave function, to the 
same value E = —7.17 ± 0.02 K. This value, in spite of 
the small size of the considered system, is close to the 
experimental^! result E = —7.14 K. SWF converges very 
quickly, in fact t = 0.05 K _1 is already enough to ensure 
convergence. CWF instead, requires a three times larger 
imaginary time, i.e. r = 0.15 K . This was somehow 
expected, since SWF is presently the best available vari- 
ational wave function for 4 HeJ^ Nevertheless, the quick 
convergence of also CWF is a really remarkable result. 
In fact, this means that PIGS efficiently includes the ex- 
act interparticle correlations through the imaginary time 
projections, without any need of importance sampling. 
Then, the choice of a good wave function, within the 
PIGS method, becomes a matter of convenience rather 
than of principle, since better initial wave functions only 
allow for a smaller total projection time r, and thus less 
CPU consuming simulations. 

This convergence is confirmed also by the radial distri- 
bution function g(r) and the static structure factor S(k). 
For such quantities, the convergence rate is found to be 
similar to the energy one. In Fig. [3] we report the radial 
distribution function g(r) obtained by projecting both a 
SWF and a CWF at different imaginary time values. For 
t > 0.05 K , SWF results at different r are indistin- 
guishable within the statistical uncertainty (see Fig. [3^,). 
In fact, with SWF the exact result is reached within very 
few projection steps and then it is no more affected by 
further projections. As already pointed out, also CWF 
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Figure 3: Radial distribution function g(r) for bulk liquid 
4 He computed in a cubic box with N = 64 at the density 
p = 0.0218 A" 3 with the PIGS method, a) g(r) obtained by 
projecting a SWF for r = 0.00, 0.05 and 0.25 K" 1 . The r = 
0.00 result corresponds to the variational SWF estimate of 
g{r). b) g(r) obtained by projecting a SWF for r = 0.25 K" 1 
and a CWF for r = 0.25 K _1 . In the inset a zoom of the first 
maximum region, c) Ag r (r) = <?swF( r )~9cwp( r ) at different 
r values, where jj W (r) is the g(r) computed by projecting 
a SWF for an imaginary time equal to r, and jJ WF (r) is the 
same but by projecting a CWF. Note the smaller scale on the 
vertical axis 



displays a fast convergence, as shown in Fig. (3J:, where 
&g T {r) = SswfM ~ ScwfM is shown. For increasing 
r, Ag T evolves toward a flat function, meaning that the 
systems described starting from the two different wave 
functions, i.e the strongly correlated quantum liquid of 
SWF and the ideal gas of CWF, are evolving into the 
same quantum liquid, which is the best reachable rep- 
resentation of the exact ground state of the simulated 
system. The same conclusion is inferred from the evolu- 
tion of the static structure factor S(k), which is plotted 
in Fig. |H 



2. PIGS results from a wrong initial function 

In order to put a more stringent check on the PIGS 
method ability to converge to the exact ground state 
without any variational bias, we have considered also a 
wrong initial wave function by projecting a GWF. Thus 
at the starting point of the imaginary time path there is 
now a strongly localized Einstein crystal. We find, even 
in this case, that the energy converges to the same value 
as before (see Fig. [5]). Thus PIGS is able not only to 
drop from the initial wave function the wrong informa- 
tion of localization, but also to generate at the same time 
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Figure 4: Static structure factor S(k) for bulk liquid 4 He 
computed in a cubic box with N = 64 at the density p — 
0.0218 A" 3 with the PIGS method, a) S(k) obtained by 
projecting a SWF and a CWF for r = 0.05 K" 1 . b) S{k) 
obtained by projecting a SWF and a CWF for r = 0.40 K _1 . 
c) AS T (k) = Ss WF (k) — Scwp(fc) at different r values, where 
S'g WP (fc) is the S(k) computed by projecting a SWF for an 
imaginary time equal to r, and Sc WP (fc) is the same but by 
projecting a CWF. Note the smaller scale on the vertical axis. 



the correct correlations among the particles. GWF needs 
r = 0.5 K" 1 to converge, which is ten times larger than 
the SWF value. 

Again this convergence is confirmed also by the radial 
distribution function g(r) and the static structure factor 
S(k). In Fig. Owe report the radial distribution function 
g(r) obtained by projecting a GWF at different imag- 
inary time values compared with the ones coming from 
the projection of SWF. It is evident that small imaginary 
time is not enough to leave out the wrong information in 
the GWF. For lower r values, there are still reminiscences 
of the starting harmonic solid, which are progressively 
lost as the projection time increases. This is made clearer 
in Fig. [SJd where we plot the difference Ag T (r) , at fixed 
imaginary time r, between the g(r) computed by pro- 
jecting the SWF and the one obtained by projecting the 
GWF. A similar behavior is observed in the evolution 
static structure factor S(k), plotted in Fig. [5] For the 
GWF, the Bragg peak shown at small r values (Fig. [6^,), 
which is typical of the solid phase, becomes lower and 
lower as the projection time is increased (Fig. [BJa), until 
convergence is reached (see Fig. [BJ:). 

From the plot of the energy per particle vs. the to- 
tal imaginary time r it is possible to estimate the over- 
lap per particle of the initial wave function on the ex- 
act ground stated By using the results in Fig. [2] we 
find that the overlap of SWF is about 99%, while the 
GWF one is about 10%. That SWF has an high overlap 
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Figure 5: Radial distribution function g(r) for bulk liquid 
4 He computed in a cubic box with iV = 64 at the density 
p = 0.0218 A" 3 with the PIGS method, a) g(r) obtained 
by projecting a SWF for r = 0.40 K" 1 and a GWF for r = 
0.50 K _1 . In the inset a zoom of the first maximum region, 
b) Ag T (r) = SswfM — AgwfM a t different r values, where 
PswfM i s the g(r) computed by projecting a SWF for an 
imaginary time equal to r, and jQ W (r) is the same but by 
projecting a GWF. Note the smaller scale on the vertical axis. 
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Figure 6: Static structure factor S(k) for bulk liquid 4 He 
computed in a cubic box with N — 64 at the density 
p = 0.0218 A" 3 with the PIGS method, a) S(k) obtained 
by projecting a SWF and a GWF for r = 0.05 K" 1 . It is 
evident in the GWF result the presence of the Bragg peak. 
Note the logarithmic scale, b) S(k) obtained by projecting 
a SWF for r = 0.40 K" 1 and a GWF for r = 0.50 K" 1 . 
The Bragg peak is no more present in the GWF result, c) 
AS T (k) = Sg WF (fc) — SgwfM at different r values, where 
SswfW is the S(k) computed by projecting a SWF for an 
imaginary time equal to r, and S , Q WF (fc) is the same but by 
projecting a GWF. Note the change of the vertical scale. Er- 
ror bars are smaller than the used symbols. 
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Figure 7: One-body density matrix pi obtained from PIGS 
simulations for liquid 4 He at the equilibrium density p — 
0.0218 A" 3 by projecting a SWF, a JWF and a GWF for 
an imaginary time r = 0.30, 0.40 and 0.80 K _1 respectively. 
The dotted line indicates the condensate value no = 0.069 
obtained from an independent PIGS simulation. 26 

with the ground state is not a surprise; it was qualita- 
tively expected since SWF is presently the best available 
wave function for 4 HeJ^ However a 99% overlap is re- 
ally remarkable and provides a further argument on the 
goodness of SWF. On the other hand, a poor overlap 
of GWF was somehow expected, since the parameter C 
was chosen to strongly localize the atoms of the bulk liq- 
uid around fictitious equilibrium positions on a regular 
lattice. 



Dirac delta function, which indicates a macroscopic oc- 
cupation of a single momentum state, i.e. Bose-Einstein 
condensation. 

The exact p\ can be obtained in PIGS simulation by 
substituting ipo in (|17[) with ip T with r large enough. This 
corresponds to the simulation of a system of N — 1 lin- 
ear polymers plus a polymer which is cut into two halfs, 
called half-polymers, one departing from r and the other 
from r* '. Thus p\ is obtained by collecting the relative dis- 
tances among the cut ends of the two half-polymers dur- 
ing the Monte Carlo sampling. The present computation 
of pi has been obtained by implementing a zero temper- 
ature version of the worm algorithm^ We have worked 
with a fixed number of particles and not in the grand 
canonical ensemble, similarly to what has been done at 
finite temperature in Ref. 23. In practice this corresponds 
to a usual PIGS calculation of p\ where "open" and 
"close" moves have been implemented 8 in order to visit 
diagonal and off-diagonal sectors within the same simula- 
tion. The advantage of doing this does not come from the 
efficiency of the worm algorithm to explore off-diagonal 
configurations, because similar efficiency is obtained with 
PIGS when "swap" moves are implemented. 11 - The ben- 
efit in using a worm-like algorithm here instead comes 
from the automatic normalization of p\ which is a pecu- 
liarity of this method^ In Fig. [7] we report p\ obtained in 
PIGS simulations of bulk liquid 4 He at p = 0.0218 A~ 3 
by projecting either a SWF, a JWF and a GWF. All the 
simulations give the same result, shown in Fig. [7] which 
turns out to be compatible with the recent estimate ob- 
tained with PIGS given in Ref. H of n = 0.069 ± 0.005. 



B. Solid 



3. Off-diagonal properties 

Besides the diagonal ones, also off-diagonal proper- 
ties, such as the one-body density matrix, are acces- 
sible within PIGS simulations. The one-body density 
matrix pi(r,r l ) represents the probability amplitude of 
destroying a particle in r and creating one in f 1 . Its 
Fourier transformation represents the momentum distri- 
bution. In first quantization p\ is given by the overlap 
between the normalized many-body ground state wave 
functions ^o(-R) and iJ}q(R'), where the configuration 
R' = (f 1 , r*2, ■ ■ ■ j rjsr) differs from R — (r, Ti, . . . , r/v) only 
by the position of one of the N atoms in the system. If 
ipo(R) is translationally invariant, p\ only depends on the 
difference \?— r*|, thus 

Pl (r-f*)=N J df 2 ...dr N %{R)MR')- (17) 

The Bosc-Einstcin condensate fraction no is equal to the 
large distance limit of pi(f — r*). In fact, if p\ has 
a nonzero plateau at large distance, the so called off- 
diagonal long-range order (ODLRO), its FT contains a 



We have performed the computation at density p = 
0.0313 A~ 3 , where 4 He is in the solid phase, by projecting 
a SWF and a CWF. Our results for the energy per parti- 
cle are plotted in Fig. [5] as a function of r. In both cases 
we find convergence to the value E = —5.34 ± 0.02 K. 
Even in this phase the convergence of SWF is faster, 
being r = 0.05 K _1 enough to reach convergence. In 
the case of CWF convergence is reached only for a much 
larger imaginary time r = 0.80 K _1 . 

Also in this case convergence is obtained for the radial 
distribution function and for the static structure factor, 
reported in Fig. [9] and Fig. [10] respectively. From Fig[9^ 
it is evident that SWF has reached the true ground state 
with few projection steps, since the results for g(r) at 
r = 0.05 K _1 and t = 0.80 K _1 are indistinguishable. 
The evolution toward the correct ground state of the pro- 
jected CWF is instead detectable. The presence of the 
crystalline structure is mainly evident in the static struc- 
ture factor, where a Bragg peak grows with increasing t 
(see Fig.HOk.b). The emerging of the correct solid struc- 
ture by projecting a really poor wave function such as the 
CWF is made evident by the trend toward a flat function 
of the differences Ag T (r) and AS T (k) plotted in FigJSJi 
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Figure 8: Energy per particle £ as a function of the total 
projection time r obtained from PIGS simulations of an fee 
4 He crystal at the density p = 0.0313 A~ 3 by projecting a 
SWF (filled circles) and a CWF (open circles). Dashed line 
indicates the convergence value E — —5.34 ± 0.02 K. 



estimator for the Hamiltonian (Eq. (JSJ with O = H) 
implying the necessity of an analytical formulation fo r 




and FigfTUb respectively. 

IV. CONCLUSION 

In this work we have studied with the Path Integral 
Ground State method diagonal and off-diagonal proper- 
ties of a strongly interacting quantum Bose system like 
the bulk liquid and solid phases of 4 He. We have obtained 
convergence to the ground state values of quantities like 
the total energy, the radial distribution function, the 
static structure factor and the one-body density matrix 
projecting radically different wave functions: equivalent 
expectation values in the liquid phase have been obtained 
using as initial wave function a shadow wave function, a 
Gaussian wave function with strongly localized particles 
of an Einstein solid without interparticle correlations and 
also a constant wave function where all configurations of 
the particles are equally probable. Similarly in the solid 
phase equivalent expectation values have been obtained 
by considering a shadow wave function, which describes 
a solid, and a constant wave function which describes an 
ideal Bose gas. The present analysis demonstrates the 
absence of any variational bias in PIGS; a method that 
can be thus considered as unbiased as the finite temper- 
ature PIMC. This remarkable property comes from the 
accurate imaginary time propagators, exactly the same 
used with PIMC, that do not depend on the initial trial 
state. It remains true that the use of a good variational 
initial wave function greatly improves the rate of con- 
vergence to the exact results. Moreover, very poor wave 
functions have also the drawback of requiring the direct 



Figure 9: Radial distribution function g(r) for bulk solid 
4 He computed in a cubic box with N = 32 at the density 
p = 0.0313 A" 3 with the PIGS method, a) g(r) obtained by 
projecting a SWF for r = 0.05 and 0.80 K _1 . b) g(r) ob- 
tained by projecting a SWF and a CWF for r = 0.80 K _1 . In 
the inset a zoom of the first maximum region, c) Ag T (r) = 
js W (r) — jc WF (r) at different r values, where jg WF (i l ) is 
the g(r) computed by projecting a SWF for an imaginary 
time equal to r, and jJ WF (r) is the same but by projecting a 
CWF. 



the small imaginary time propagator G or an accurate 
knowledge of its derivatives. 

We have addressed here only the case of a realistic in- 
teraction potential among Helium atoms. However one 
can reasonably expect that this conclusion holds even for 
very different kinds of interaction, once an accurate ap- 
proximation for the imaginary time propagator is known 
(for example hard-spheres 27 or hydrogen plasma 28 ). As 
far as pathological potentials like the attractive Coulomb 
one are concerned, PIGS would suffer the same limita- 
tions of PIMC if inaccurate approximations of the prop- 
agator were used^ 
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Figure 10: Static structure factor S(k) for bulk solid 4 He 
computed in a cubic box with TV = 32 at the density p — 
0.0313 A -3 with the PIGS method, a) S(k) obtained by pro- 
jecting a SWF and a CWF for r = 0.05 K" 1 . b) S{k) obtained 
by projecting a SWF and a CWF for r = 0.80 K _1 . The black 
dots are under the red ones, c) AS T (k) = Sswf (k) — Sqwf CO 
at different r values, where Sg WF (fe) is the S(k) computed 
by projecting a SWF for an imaginary time equal to r, and 
Sc WF (fc) is the same but by projecting a CWF. Error bars are 
smaller than the used symbols. Notice the logarithmic scale 
in panels a) and b). 
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